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Abstract 


The baste governing equation* for the second-order three-dimensional hy- 
personic thermal and chemical nonequilibrium boundary layer are derived 
by means of an order-of-magnitnde analysis. A two-temperTure IceTis 
mplemented into the system of boundary-layer equations by simpSthe 
ather complicated general three-temperature thermal gas model P The equa- 
tions are written in a surface-oriented non-orthogonal curvilinear coordinate 
system, where two curvilinear coordinates are non-orthogonal and a third 

imum 15 n ° rmaI to the surface - The equations are described with min- 
~ ° f tCnSOr eXpTe f 3ns arising from the coordinate transformation 

will be suThTT*^ 1 T ^ ieaderS - The SCt ^uations obtained 
wiU be suitable for the development of a three-dimensional nonequilibrium 

boundary-layer code. Such a code could be used to determine economicX 
he aerodynamic/aerothermodynamic loads to the surfaces of hypersonic ve- 
hicles with general configurations. In addition, the basic equations for three 
dimensional stagnation flow, of which solution is required as an initial value 

al ' ntegration of the boundary-layer equations, are given 
along with the boundary conditions, the boundary-layer parameters and the 

rioTT 7:; m ; tCU ” g pr °“ d " r '- for the chemfoJreat 

noneouffibri^ ^ ltermod >'" amic transport properties in the thermal 
nonequilibrium environment are explicitly given. 
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Nomenclature 


A 

= 

vector quantity 

a 

= 

determinant of surface metric tensor 

£ 0 , a° 

= 

covariant and contravariant surface base vectors 

®a/3 1 ® ^ 

= 

surface metric tensor components 

ba0 

=r 

component of covariant surface curvature tensor 

c. 

= 

mass fraction of species s 

C P 

= 

frozen specific heat of mixture at constant pressure 

C T*,' 


specific heat of species s at constant pressure 
for translational mode 

C pel 

= 

specific heat of mixture at constant pressure 
for electronic excitation 

c pel,$ 

= 

specific heat of species s at constant pressure 
for electronic excitation 

C P r,» 

— 

specific heat of species s at constant pressure 
for rotational excitation 

Cpv 

= 

specific heat of mixture of vibrationally excited molecules 


= 

specific heat of molecular species s at constant 
pressure for vibrational excitation 

Cvtl t 9 


specific heat of species s at constant volume 
for electronic excitation 

CvV,9 

= 

specific heat of species s at constant volume 
for vibrational excitation 

D. 

= 

effective diffusion coefficient of species s 

D„ 

=r 

binary diffusion coefficient for s-r species 

d : 


ambipolar diffusion coefficient of ionic species s 


“ 

average vibrational energy per unit mass of molecular species s, 
which is created or destroyed at rate w s 

e 


mixture energy per unit mass 


~ 

mixture electronic energy per unit mass 

&el,i 

— 

electronic energy per unit mass of species s 

e. 

= 

energy per unit mass of species s 

e°. 

= 

energy of formation of species s 


= 

mixture vibrational energy per unit mass 
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vibrational energy per unit mass of species s 
equilibrium vibrational energy per unit mass of 
species s at translational temperature T 
first ionization energy of species s per mole 
scalar quantity 

determinant of general metric tensor 

covariant and contravariant general base' vectors 

covariant and contravariant general metric tensor components 

mixture enthalpy per unit mass 

mixture electronic enthalpy per unit mass 

electronic enthalpy per unit mass of species s 

enthalpy per unit mass of species s 

heat of formation of species s 

mixture vibrational enthalpy per unit mass 

vibrational enthalpy per unit mass of species s 

Jacobian of transformation 

Boltzmann constant 

backward reaction rate coefficient for reaction p 
forward reaction rate coefficient for reaction p 
equilibrium constant for reaction p 
body length 

Lewis number of species s 

Lewis number of vibrationally excited molecular species s 

mixture mass 

mass of species $ 

molecular weight of species s 

freestream Mach number, or mixture molecular weight at freestream 

molar rate of production of ionized species s 

per unit volume by electron-impact ionization 

pressure 

Prandtl number 

Prandtl number for mixture of vibrationally excited molecules 

heat-flux vector 

wall heat-transfer rate 

ordinary gas constant 
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ft 

Re 

S. 

T 

T a 

T th 

T v 

T Vi ,h 

Tij , 1* 
Tr tf 
l 7cc 

U, V, w 

U, V, W 
w, 

V. 

v' 

V, 

y, z 

x { 

y v 

y» 

a*, 

0f 

7 . 

It 

r J 

1 ki 

8 

8 a 

3 

Ag\ AW 


— universal gas constant 
= Reynolds number 

= bridging function in vibrational rate equation 
= translational-rotational temperature of heavy particle 
= geometrically averaged temperature for dissociative reactions 
post-shock translational-rotational temperature 
vibrational-electronic-electron temperature 
= post-shock vibrational temperature 
= covariant and contravaxiant stress tensor component 
= reference temperature 
= freestream velocity 

— Cartesian velocity components 

contravaxiant velocity components in a curvilinear system 
mass production rate of species s per unit volume 
= velocity vector 
= ith component of velocity 
= diffusion velocity vector of species s 
= Cartesian coordinates 
= general non-orthogonal coordinate 
~ Cartesian reference coordinate 

— molar fraction of species s 

= component of contravaxiant base vector 
= component of co variant base vector 
= molar concentration of species s per unit mass 
= total molar concentration of mixture 
= Christoffel symbols of the second kind 
= boundary-layer thickness 
= mass flow displacement thickness 
= Kronecker tensor 

modified collision integrals for the collision pair s — t 
characteristic energy for diffusion model for molecular species s 
= frozen thermal conductivity for translational-rotational 
energy of heavy particles 
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frozen thermal conductivity of electrons owing to 

collisions between electrons and all particles 

frozen thermal conductivity of electrons owing to 

collisions among electrons 

thermal conductivity for rotational energy 

frozen thermal conductivity for translational energy 

thermal conductivity of vibrational energy 

owing to collisions between molecules and all particles 

thermal conductivity of vibrational energy 

owing to collisions among molecules 

second coefficient of viscosity 

viscosity coefficient of mixture 

reference viscosity 

curvilinear surface-oriented coordinates 
mixture mass density 
mass density of electrons 
overall mass density of molecules 
mass density of species s 
limiting collision cross-section 

collision limited vibration al relaxation time for molecnlar species s 
combined vibrational relaxation time of molecular species s 
average vibrational relaxation time for species s 

vibrational relaxation time for species s with collision partner r 
dissipation function 

transformed normal coordinate, (= (/6) 



Subscripts 
e = electron 

ed = boundary- layer outer edge 

el = electronic excitation mode 

* — index of the nodes in x * direction 

j = index of the nodes in direction 

k = index of the nodes in x * direction 

r = species r, or rotational mode 

ref = dimensional reference quantity 

s = species 5 , or stagnation point; 

1 = TV, 2 = O, 3 = N 2 , 4 = 0 2 , 5 = NO, 6 = N+, 

7 = 0 + , 8 = JV a +, 9 = 0+, 10 = NO + , 11 = e" 
sh = post-shock condition 

v = vibrational mode 

w = quantity at the wall 

oo = freestream quantity 

Superscripts 

— — dimensional quantity 

* = physical quantity 

Conventions 

Latin indices in tensor relationships vary between 1 and 3, Greek indices 
between 1 and 2. 
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1* Introduction 


th.7 i !V he J raP ^ adVen ‘ h “ dw " e “ d numerical algorithms in 

the last decade it is now possible to solve the Navier-Stokes (N-S) equations 

^ “°“ <l ^ b ° dy 

Jl:rr l l u ffiCUlt ’ !f n0t impossible ’ to ™P™duce experimentally 
omplex flow fields around present and future space vehicles, such as the 

Space Shuttle Orbiter, assured crew return vehicle (ACRV), and aeroassisted 
space transfer vehicle (ASTV), in a ground-based LulatL fac“hu 6 
i ecomes lg y desirable to rely on computational-fluid-dynamics (CFD) 
or hypersonic flow simulation. However, systematic application of the N-S 
solvers to the design of the space vehicles is still limited by the extremely high 
costs of numerical operation. This is especially true for solving the three- 
dimensional flow field around a complete geometry, which needs prohibitive 
amounts of computer time and storage. 

r,JI U T" CiI ‘ 00ls f° r d'rigu environments requite fast solution of many 
rod j' eseiltl ? availab| e three-dimensional methods such as a streamline 
code and an asymmetric analog code [3, 4) are too approbate and too 
restnctive for general use. Furthermore, these codes can not account for real 
g effects or surface effects which are important aspects of the hypersonic 
flight regunes. Consequently, there is a need to develop a computational code 
which is economical to use and has the required accuracy as an engineering 

J**” .°°: " a ttmctive candidate for meeting this requirement is the 

viscous-, nviscd two layer approach. By matching the boundary-layer Sow to 
the external inviscid flow about three-dimensional configurations, an entire 
effort ^ h' prcdicted witl “ reduction of the computational 

Until the mid 70’s, study of three-dimensional compressible boundary- 
ayer theory had been actively conducted [5-8] in the U.S. in relation to 
the development of hypersonic aircraft. Although these three-dimensional 
studies were extensive and involved laminar, transitional, and turbulent 
boundary-layer flows, all of them were based on classical, first-order boundary 
ayer heory. As a result, practical application was limited to smooth bod- 
le s hke conical geometries and thin wing configurations. In addition, most 
methods utilized an orthogonal curvilinear coordinate system which further 
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hindered their use for broad application. Meanwhile, research efforts then 
were shifted to the development of numerical techniques for solving the full 
Navier-Stokes (FNS), parabolized Navier-Stokes (PNS) and viscous-shock- 
layer (VSL) equations, leading to the dwindling of boundary-layer approach. 

Classical, first-order boundary- layer theory is based on the assumption of 
very small curvature compared with the boundary-layer thickness. This is 
not the case for a general three-dimensional geometry. Furthermore, classical 
boundary-layer theory can not handle second-order effects such as displace- 
ment effects which are provoked by the viscous region, entropy and total 
enthalpy gradient effects, which originate in the inviscid region of the curved 
bow shock present in supersonic/hypersonic flow, and body curvature ef- 
fects. The shock wave curvature in the nose region has a strong influence 
on downstream aerodynamic heating and on the boundary-layer edge prop- 
erties, which in turn determine the convective heating load of the vehicle [1]. 
The limitations of the first-order boundary-layer theory can be overcome to 
a certain extent by considering these second-order effects properly. 

Second-order boundary-layer theory was first studied by Van Dyke [fi- 
ll], who distinguished the second-order effects. He expanded flow variables 
asymptotically in powers of the perturbation parameter Re -1 / 2 , where Re is 
a characteristic Reynolds number. This was done separately for the outer 
inviscid and inner viscous flows. The results are then substituted into the 
Navier-Stokes equations to obtain the first- and second-order set of equations 
for each flow region. To obtain the matching conditions where the regions 
overlap, the variables of the outer flow were represented by an appropriate ex- 
pansion of its wall values. Van Dyke applied this approach to two-dimensional 
and axisymmetric flows. Recently, Aupoix et al. [12, 13] have modified the 
asymptotic expansion approach in such a way that a smooth matching of 
the inner and outer solutions is ensured. Their so called ‘defect’ equations 
include Van Dyke’s equations plus terms of higher order. An attempt to 
extend this approach to three-dimensional flow is under way. 

An alternate way to obtain the second-order boundary-layer equations 
is the order-of-magnitude analysis which has been widely used for deriving 
the simplified form of the Navier-Stokes equations. In this analysis, the N-S 
equations are first expressed in dimensionless form by normalizing all flow 
variables with reference values while all linear dimensions are normalized by 
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a characteristic length. Each term of the N-S equations is then estimated and 
those smaller than a prescribed value are dropped from the equation system. 
In the first-order theory, only terms of order unity axe retained leading to the 
well known classical boundary layer equations. In the second-order theory 
terms up to order Re -1 / 2 , which is the same as the perturbation parame-’ 
ter defined by Van Dyke [9-11], are retained. This approach is much easier 
to understand than the perturbation approach, since it follows the orderly 
simplification of the full Navier-Stokes equations thus allowing the compari- 
son with PNS and VSL equations depending on the level of approximation. 
Davis and Flugge-Lotz [14] used this order-of magnitude method comprehen- 
sively to develop the VSL equations which are the composite equations for 
the Euler equations and the second-order viscous equations. 

It must be noted that the surface curvature effect is customarily distin- 
guished from the other second-order effects such as the displacement effect 
the external entropy gradient effect and the external total enthalpy gradient’ 
All these second order effects only affect the matching of the boundary-layer 
flow with the external flow and thus influence the boundary conditions of 
either the outer inviscid region or the inner viscous region without directly 
intervening in the boundary layer equations. The displacement effect can 
be dealt with by means of either the displacement surface or the equivalent 
inviscid source distribution. The entropy and total enthalpy gradient effects, 
which originate in the inviscid flow region, are taken into account through 
the iterative viscous-inviscid matching process. 

On the other hand, the curvature effect is caused exclusively by the shape 
of the body surface. When the local boundary-layer thickness is not small 
compared with the smallest radius of curvature of the body surface, which is 
the case for any practical body geometry, the curvature produces a non-zero 
pressure gradient across the boundary layer, expressed by the body-normal 
momentum equation. Moreover, the curvature contributes other terms to 
the second-order equations and affects the boundary conditions. 

In order to predict correctly the flow field about general three-dimensional 
geometries based on the viscous-inviscid two-layer model, one must have an 
elaborate second-order boundary-layer solution, an accurate inviscid Euler 
solution and a rigorous inner-outer matching technique. Using this strategy, 
Monnoyer [15-17] has calculated the hypersonic three-dimensional laminar 
flow on general body configurations. His method is a direct extension of 
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the comprehensive European efforts [18, 19] in the development of higher- 
order boundary-layer theory. Monnoyer [15] closely followed the order-of- 
magnitude analysis of Robert [18] , who derived the second-order boundary- 
layer equations from the N-S equations in a curvilinear non-orthogonal coor- 
dinate system [20], and modified the evaluation of surface curvature. Mon- 
noyer applied the second-order boundary-layer equations for solving the three- 
dimensional flow around an ellipsoid at incidence [16], and then, flow past the 
windward side of the European space shuttle HERMES [17]. Although Mon- 
noyer has dealt with three-dimensional boundary-layer flow quite extensively, 
only primitive gas models (perfect and equilibrium gas) have been applied 
so far. Monnoyer’s work and other recent papers [21, 22] demonstrate that 
the inviscid/boundary layer approach can provide flow predictions for cases 
without separation that are nearly as accurate as Navier-Stokes solutions and 
are significantly cheaper. 

Apart from the European work, an attempt was made to develop the 
higher-order boundary-layer equations at the NASA Ames Research Center 
in the late 80 ’s. Panaras [23] derived the second-order boundary-layer equa- 
tions for three-dimensional compressible flows. The equations were written 
in a generalized curvilinear coordinate system, and were developed from the 
Navier-Stokes equations written in tensor notation by means of an order-of- 
magnitude analysis. The resulting equations are nearly identical with those 
developed by Monnoyer [15, 16]. Subsequently, Steger et al. [24, 25] mod- 
ified Panaras’s analysis to describe the boundary-layer equations in general 
body-fitted curvilinear coordinates while retaining the original Cartesian de- 
pendent variables. This alternate form of the boundary-layer equations offers 
several advantages in terms of numerical stability, and the software (gridding, 
boundary condition routines, etc) developed for many Navier-Stokes schemes 
can be readily applied. This attractive new formulation of the boundary-layer 
equations is, however, rather arbitrarily simplified when applied to sample 
calculations because of concerns about numerical stability. Very limited cal- 
culations for a perfect gas have been undertaken so far. Therefore, the general 
use of this approach remains to be seen. 

In the mid 70’s, Gershbein [26] in the former U.S.S.R. developed the sys- 
tem of equations for a three-dimensional chemically reacting laminar bound- 
ary layer in curvilinear non-orthogonal coordinate system. His set of equa- 
tions was a direct extension of the three-dimensional boundary-layer equa- 
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tions derived by Shevelev [27] for a perfect gas, which was based on first-order 
boundary layer theory. Due to the first-order analysis, Gershbein dealt with 
only flows past simple smooth blunt bodies such as an ellipsoid and a cylinder 
with a primitive gas chemistry. 

In the hypersonic flight regime, the thermal and chemical characteristics 
of air in the shock layer of space vehicles are altered in ways that affect 
the thermodynamic and transport properties, the chemical reaction rates, 
and the radiation properties. These properties, in turn, influence the char- 
acteristics of the flow, the shear stress, and the heat flux. For the next 
generation of aerospace vehicles, which usually will fly through the upper 
atmosphere of Earth, nonequilibrium aerothermochemistry and finite-rate 
surface catalysis will play an important role in determining the aerody- 
namic/aerothermodynamic loads to the surface [28-32]. 

Highly sophisticated high-temperature thermochemical gas models have 
been developed in recent years at the NASA Ames Research Center [33]. 
These real gas models have been evolved for the purpose of maturing the 
enabling technology necessary for the design of proposed aerobraking and 
aeromaneuvering space vehicle configurations such as an aeroassisted orbital 
transfer vehicle (AOTV) and transatmospheric vehicle (TAV). The remain- 
ing task, therefore, is how to simplify, without significant compromise, these 
rather complicated thermochemical gas models and to implement them into 
the engineering flow field solver such as the boundary-layer method. 

Since the type of the boundary-layer equations is a parabolic, they can be 
solved using space-marching integration methods. The Navier-Stokes equa- 
tions, on the other hand, have to be integrated with time-marching methods 
which require enormous amounts of computational time and storage. Numer- 
ical schemes for the boundary-layer equations are relatively well-established 
and tested as compared to methods for the other viscous flow equations [34]. 

The purpose of the present study is to derive the second-order hypersonic 
three-dimensional thermochemical nonequilibrium boundary-layer equations 
in generalized curvilinear coordinates. The resulting set of equations will be 
suitable for the development of a boundary-layer code which allows simulta- 
neous handling of both general three-dimensional configurations and general 
thermochemically nonequilibrium gas models in an engineering context. Al- 
though the future options include wall catalysis, ablation, radiative heat 
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transfer and turbulent flow, major focus is placed on the three-dimensional 
laminar boundary-layer flow field around practical geometries with reed-gas 
effects. The strong-viscous interactions, such as flow separation or shock 
boundary-layer interaction, will not be treated here since these problems are 
beyond the scope of the boundary-layer approach. 
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2. Governing Equations 


2.1 Coordinate System 

The selection of a coordinate system to describe the boundary-laver eo ,,». 

ca^y^ecttfi”* -dimensional body is of major importance and can signifi- 

In i k , “ d a PP Ucab ^y of the final computer cod! [7] 

In addition, the boundary-layer approximation is impractical to implement 

“ s y«tem other than a surface-oriented coordinate system sZ* 
an orthogonal system may cause a number of 

lengthy interpolation procedures [6]. Consequently, the natural choice for th! 
connate system of general three-dimensional boundary-layer equations is 
a surface-oriented non-orthogonal curvilinear coordinate system wh^ two 
curvilinear coordinates are locally parallel to the surface and a third coordi 
nate is normal to the surface in order to implement a thin boundary-layer t 
suction. In this coordinate system, surface-normal derivative terms whkh 
Play^an important role in the boundary-layer equations, can be eLTly r^ 

W ThC SUrface -° rient f d non-orthogonal curvilinear coordinate system has 
en widely used in the works of European ktouds fio.onl i , , 

named it a surface-oriented locally monoclimc coordiiate sjLtem.'rW utd 
extens.ve tensor eoncepts, which are compHcated enough to scare off the un 
farmhar reader to derive the boundary-layer equations and to determtae Z 
metnc propert.es of the body surface. Furthermore, their unique u ””ndex 
notatron and so-called ‘shifters’ may cause unnecessary confusion. On the 

the? t a f d ’ P J I ? ra ® described tbe system of equations in such a way that 
he interface between the analytical partial-differentiation notation Jd the 

^TdThe H Y CaSily “ nderstandable * Therefore, the tensor formulation 
mid the order-of-magnitude analysis of Panaras [23] are used in the present 

study along with the boundary condition treatment, boundary-layer parame- 
ers and inner-outer layer matching procedure of the European group^-igl 

^Tw^ 00 ’ fT ever> the usual notation of fluid d y namics is employed 
with minimal use of tensor expressions. ^ 


2.1.1 Base Vectors and Metric Coefficients 
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The general non-orthogonal, curvilinear surface-oriented coordinate sys- 
tem a :’(£,f7,C) i s depicted in Fig. 1 where the Cartesian coordinate system 
y'\x,y , z) is used as a reference. The coordinate system consists of two 
curvilinear coordinates £ and tj locally parallel to the surface and a third 
coordinate ( locally orthogonal to the other two coordinates and rectilinear. 
It must be noted that the complete space above the surface can be uniquely 
described by these coordinates if the surface is convex but only partly if the 
surface is concave [15]. 


In the present paper, only the essential geometric relations, which result 
from coordinate transformation and which are presented in tensor notation, 
are described. The detailed tensor analysis applied to the basic equations of 
fluid dynamics can be found in [15, 19, 20]. In tensor notation, superscripted 
indices, e.g. it*, are called contravariant tensors and subscripted indices, e.g. 
Ui, covariant tensors. In Fig. 1, e,< is an unit vector in the Cartesian reference 
coordinate system, is a covariant vector at the surface (£ = 0) and g. is a 
general covariant base vector (oflf-surface points). The general covariant base 
vector g. in this system is defined by the transformation [20] 


gi 


dy v 

dx i§i ' 


& 


= 1,2,3) 


( 1 ) 


where the notation /3f is defined by 


./ dy * 

Pi = 


dx x 


( 2 ) 


In the present notation, the Einstein summation convention is used through- 
out, e.g.: 

PI &i> = P\ &i> + Pi £ 2 ' + Pi £ 3 ' 


In a general non-orthogonal curvilinear coordinate system, the base vector g. 
is neither perpendicular nor of unit length contrary to the case for a Cartesian 
coordinate system. 

Similarly, the covariant base vector at the surface is defined by 


a, = /?*£,; (*,*' = 1,2,3) 


( 3 ) 


where the superscript *•’ denotes a quantity at a point on the body surface. 
The contravariant base vector g x is defined by 


i dx< i 

£ = = a *'^ 


(i,i'= 1,2,3) 


( 4 ) 
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where aj/ is defined by 


a’,= 


ax' 

<V 


The contravariant base vector must satisfy the relation 
where 6 * is the Kronecker tensor 


( 5 ) 


( 6 ) 


<? = 


1 0 0 
0 1 o 
0 0 1 


td^T “° rdi, “ te C0 ' tari “ t “<• 

Ti e covariant metric tensor components g„ is defined by 

( 7 ) 

“o'n" me ‘ riC ‘ eaS01 C °»“ ' »«h- % the orthogonal- 

9ij9 a = /g\ 

^tM^re" 6 ^ ‘° ““ ° f * h ' ™* -trie tensor 

It. 

( 9 ) 


& g ( 9jm9kn ~ 9jn9km ) 

where indices and (Z.m.n) are cyclic. 

Here, g is the determinant of the metric tensor 


9 = l9ijl = 


9u 9i2 9i3 
921 922 g 23 
931 932 933 


1 


Therefore, the Jacobian of the transformation from /(*,*,*) to *<({ |jy .f) 


s 0(1-7. C) 

d ( x ,y,z) 


6 t v 6 

Vx T)y Vz 

Cv c* 


= la- 
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is related to the determinant g by the relation 


J = 


1 

V9 


( 10 ) 


When differentiating tensors in terms of curvilinear coordinates, certain com- 
binations of partial derivatives of metric tensor components (the Christoffel 
symbols) play an important role. The Christoffel symbols of the second kind 
are defined by 

I* = r fc = — • g h 
dxi 2- 

rj- £ire further denoted by 


( 11 ) 


r k _ 1 Jk/ fyhj , _ ®9ij \ 

ii ~2 9 K dx i dx> dx l) 


(12) 


The Christoffel symbols may be interpreted in terms of the variation of base 
vectors with respect to the coordinates. The Christoffel symbols are zero in 
a Cartesian coordinates system. 

In the present coordinate system, the base vectors g* are space-dependent 
and g„ is an unit vector, thus gij takes the form: 

— 3 


9ij = 


9u 9n 0 

<712 922 o 

0 0 1 


and the deter minan t g of the metric tensor becomes 

9 — 5 ii 022 — ( 512) 2 ( 13 ) 

where the off-diagonal term 312 equals zero for orthogonal coordinates. For 
explicit expression of the metric tensor components, see Appendix A. 

The metric tensor at the surface is a special case of 50 [19]: 


a a0 (x a ) = gij(x a ,x z = 0); (a,/3 = 1,2) 


thus, the surface metric tensor becomes 




an ai2 

ai2 a 2 2 


( 14 ) 


ie 



The covan ant curvature tensor b Q0 of the surface is defined by [19] 

, da* . 

ba0 = d^'^ ; ( a ’/ ? = 1 > 2 ) 

With equation(3), b a0 can be expressed by 

\Qk’ 


( 15 ) 


b a0 = y - f ] 

' k w 0a "[t., 


Similar simplifications can be applied to the Christoffel symbols. 

^3a — It 3 = r 33 = Ij, = 0 (a = 1, 2) 


(16) 


Equation^ 16 ) states that any Christoffel symbol vanishes if more than one 
of its indices is a 3 (normal coordinate). For details, see Appendix B. 

2.1.2 Physical Flow Quantity 


nof , no “- Carte «“ coordinate systems, in general, tensor components do 
ot have the correct physical dimension even if they are orthogonal. This is 
ue to the fact that in such coordinate frames the base vectors are functions 
he coordinates x and are not necessarily dimensionless unit vectors as 
m a Cartesian system [20, 35]. Any vector v can be written as the sum of 
its components with respect to an arbitrary set of covariant (contravariant) 
base vectors g.(g ) referred to a particular point in space: 

TL = v% = vtf (17) 

where v x is a contravariant velocity component and v, is a covariant velocity 
component. ^ 

The contravariant physical velocity component v*‘ is defined by 


v = v*‘e. 


where e,- is an unit vector parallel to the covariant base 


(18) 


vector Since 


e. = 


y/ii m Si y/dii 


(19) 
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Then 


v*' = V& (20) 

In similar fashion, the covariant physical velocity component v* is given by 


v* = \Jg"vi 


( 21 ) 


2.1.3 Covariant Derivatives 


When differentiating a vector in terms of curvilinear coordinates, not 
only the vector itself but also the corresponding base vectors have to be 
differentiated as well since they are space- dependent. 

The partial derivative of a vector y_ consists of two parts [20]: 

dy_ 

dxi 


wherein 


is called the covariant derivative of the contravariant vector component v\ 
The covariant derivative of stress tensor component T is given by the rela- 
tion: 

T% = + Itr 1 ' + iiT* (24) 

In the case of a scalar /, the covariant derivative reduces to the partial 
derivative with respect to the coordinate 


■^ +v % 



(22) 


(23) 


/b' = 


3/ 

dxi 


(25) 


The divergence of a vector A is defined by [23] 

BA' 

divA=j- + T' ijA > 


(26) 
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Since 


( 27 ) 


(28) 

(29) 


2.2 Dimensional Form of the Navier-Stokes Equations 

In the present paper, we deal with thermal and chemical nonequilibrium 
air flow with 11 species ( N , O, N 2 , 0 2 , NO, N + , 0 + , IV 2 + , NO + , 
e~). To reduce the complexity of the general three-temperature thermal 
model developed by Lee [29], the two- temperature approach proposed by 
Park [36, 37] is adopted in the present study. For the details of general basic 
equations based on the three-temperature model for hypersonic thermochem- 
ical nonequilibrium flow in Cartesian coordinates, see Ref. 29. According 
to Park’s two-temperature model [36, 37], translational-rotational tempera- 
ture T and vibrational-electronic-electron temperature T v are necessary to 
describe the system of governing equations and to express the chemical re- 
action rates and thermodynamic and transport properties. The concept of 
a two-temperature model originates in the fact that the energy transfer rate 
between the translational mode of free electrons and the vibrational mode of 
molecular nitrogen (e-V process) is much faster than the rate between trans- 
lational and vibrational modes of nitrogen molecules (T-V process) [38, 39]. 

The necessary set of basic governing equations consists of overall mass, 
overall momentum, vibrational-electronic-electron energy, overall energy and 
species mass conservation equations plus the complementary equations of 
state, which are required to close the system of equations. In deriving the 
equations, diffusion due to the species concentration gradient is considered, 
while thermal and pressure diffusion are neglected. 
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The dimensional form of the Navier-Stokes equations for steady ther- 
mochemical nonequilibrium flow written in general curvilinear coordinates 
system x* are given by [20, 23, 29, 40] 


Continuity: 


Momentum: 


W (^v5) = ° 


dv* 


(30) 

(31) 


P^^- + pr jk ^ = T% 

where the stress tensor T' 3 and its covariant derivative T* 3 | ; - are given by 

T j = (-p + Jv k | k )g^ + « + ^In) (32) 


1 1 


;fc 


+ +<rB'u] 

+r; 1 ,^ s i, ^*u+s‘vi„) 


(33) 


Vibrational-electronic-electron energy: 


i&tv _ , dtp L. 1 d f /—/—! v 

~~ + jkV ^ + y/g dx i ^ y ^^ K ‘ v + Ke ^ 9 


e v,# e v,# j ^sh 


, — /-» V,# j J 


dx 3 

v 1 5, 


i=M 




10 


10 


+ 3 JS.WT - T.) £ £ "..A. + £ 

r=l ikZ >* i=6 


( 34 ) 


t=M 


Overall energy: 
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dh _ dp _ i a ,, ar , 
fr-SS-'SS = * + 


3F <9F 
1 5 


+ 


!v5(«. + *. V'Url + ^sK^EW^) (35) 


V /J83r ,vav " ' "■" 83P' ' ^58? 

where the dissipation function is given by 

^ = XrljU'l, + + S i ”E i |n)s J '»B*] 


(36) 


Species continuity: 


_^dc, i d ,_-=- dc,, , . 

pv dr ~ y/gdr^ pDt dx i ^ W,] - 1 > 2 »-”» 11 ) ( 37 ) 


Equation of state: 


10 W m 


Equation of state for electrons: 


P' = %pC e 


Jf' 


(38) 


(39) 


In the above equations system, an overbar l - ’ denotes a dimensional quantity. 


The heat-flux vector <f, in tensor notation, can be expressed as 


-■ _ ijdT f _^_. ij dT v “ „ ydC, 

q =-*h9 3 1 ^-(Kv + K e )g J ^-pJ2h.D.g i - 




(40) 


2.3 Normalized Form of the Boundary-Layer Equations 

The second-order boundary- layer equations are derived in a system of co- 
ordinates where the surface coordinates are denoted by (x 1 ,® 2 ) = {£,rj), and 
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the normal one by x 3 = (. In order to carry out the order-of-magnitude ana- 
lysis, the independent variables 77, £ and the flow variables are normalized 
as follows: 


</>v 

II 

H 


^►1 

II 

V 


„ r w 

U = T’ 

U OO 

w = 

JJ 

u 00 

fl __L 

n — I*- 

V- V 

p — w > 

Poo 

Pt ~ Pc’ 

P _ 7 T 2 

Poo U 00 

Pe 

II 

•> 

T 

T - , 

Poo & 00 

lv — T 

1 rtf 


II 

U OO 

h. 

11 

h v > = =r> 

^ 00 

c — is. 

C P>* F,oo’ 

fl = 

>«/ 

A - - 





h ~ iW^-’ 


Mre/^oo 

< 

/f — 

k! 

K 

Pref C poo 

fv# ^ 1 

e Pr*/ c P«’ 

K e 

Pref C poo 

w,L rt} 

)• 

1 

E* 

D v , 


to 

■J 

1 

c * a ~ n TJ > 

Poo U 00 

vl_ 

Ei,. 

X - S'*. 


Kr L "f 

mJJ 2 J 

T * “ ’ 

^ cr 

u » 


q w = 


Vr'fK 


where the overbar indicates dimensional quantities, subscript ‘ref’ refer- 
ence quantities and ‘oo’ freestream quantities. 

In the above equations, the velocity vector components U , V and W are 
the contravariant velocity components, which are connected to the Cartesian 
velocity components u, v and w by the components of the contravariant base 
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vectors 


U = ZxU + tyV + ( X W 
V = T) X U + T]yV + T) Z W 

W = ( m U + Cv v + C* W 


( 41 ) 


Furthermore, they are related to the physical contravariant velocity compo- 
nents U*, V* and W * by means of the diagonal elements of the metric tensor 
gij [equation(20)]. 

U* V* 

U = 4 =, V = 4 = (42) 




y/922 


and since 533 = 1, 


W -W* 


(43) 

On the other hand, the derivative of a scalar quantity such as the pressure 
is related to the Cartesian pressure derivatives by the chain-rule [23] 


P( = X (P* + VtPv + H Pz 


Pr, X T)Px "I" y v Pv Zf)Pz 

Pa = x cP* + y<p v + z <p* (44) 


The second-order, three-dimensional, compressible boundary-layer equa- 
tions have been derived previously from the non-dimensional Navier-Stokes 
equations by an order-of-magnitude analysis in non-orthogonal locally normal 
surface-oriented coordinate system [15-18,23]. In addition to the boundary- 
layer thickness parameter 6, which is the order of Re -1 / 2 , Monnoyer [15] 
introduced the following curvature parameter k: 

k = mox(| Ki I, | K 2 I) (45) 

where Ki and K 2 are the principal curvatures of the surface (see Appendix 
C). According to Monnoyer’s analysis, first-order theory is applied when 
O(k) < 0(1), which is equivalent to the classical boundary-layer assumption 
0(1 /k) 0(8) which states the local radius of curvature of the surface is 

much larger than the boundary-layer thickness. In this case, the variation 
of the metric coefficients across the boundary layer is negligible. On the 
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other hand, if 0 ( 1 ) < 0(k) < 0(1/8), then the variation of the metric coeffi- 
cients with distance from the surface is no longer negligible, and the second- 
order terms must be considered. This second-order condition of Monnoyer is 
equivalent to the one in Panaras’ analysis that the normal derivatives of the 
metrics axe assumed to be 0 ( 1 ) < dg^/dx 3 < 0(1/8), (i,j = 1 , 2 ) for large 
surface- curvature [23]. As mentioned in 2 . 1 , the order-of-magnitude analy- 
sis of Panaras is used here mainly because of its simplicity along with the 
methods of the European groups for boundary conditions, boundary-layer 
parameters and inner-outer layer matching [12-19]. 

The order-of-magnitude analysis of Panaras is applied to the system 
of basic governing equations (30)-(37). The resulting second-order, three- 
dimensional, laminar boundary-layer equations for a multicomponent two- 
temperature chemically reacting gas mixture in the curvilinear coordinate 
system ((,7],() are given by 


Continuity: 


^-momentum: 


SjpU^g) d(pV^g) 9{pW^g) 

9 ( ft) 8 ( 


+ pv// + P w^ + pv'rl + tyuvrl , + p v*t\, 

_ 922 9jp ffl2 9p 

9 g dr] 

1 d dU 11 dUd^g 
+ Re dC* 1 dC + AC 

1 1 d ! J ft du ! ft dV 

+ + yr23)] + Re d£ + Re «9C 


(46) 
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(47) 


77 - momentum: 

pUjr + pV ^ + pW + pU 1 It, + 2pVVT\, + pV 2 r], 
_ 912 _ 9u ?L 

~ 9 d( g dr/ 

1 d , dV^ , 1 1 JVdjg 

+ Re dC ^ d( ' + Re^/_d( d( 
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_l 1 L_ — f /nn(UT 2 4- VT 2 4- — ^~ r 2 -i- ii — T 2 

+ Re JgdC^__^ + Vr23)J + Re d£ ?31 + Re d£ 


(-momentum: 


2_ dc __ d( r* <K 
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dp 


( 48 ) 


(49) 


Vibrational-electronic-electron energy: 
pU ^L + p y^fl + = p f— 

p 3£ + 9 dr, ^ pVV dC Pe[ d( 

+ Oil + r 2 i)ir + (rj 2 + rf 2 )(Z7 + v) + (r\ 2 + r 2 2 )v] 

I I — \(k! 4- K 1 ^^ 1 4 !) ^ { K v *e) Mfgk 

+ i?e5C K v+ ej g( + -Re^C V? 3( (K 


+ 


d.^,pLe Vi , dC, 
: v2^ “51 "' v >*~a7 ’ ) + 






1 1 ^/ f pLe Vt , t ^0(7, 


+ /- e + 3 (^-w r - r.) e £ 

$—M ± ft c poo r= j 

10 

- E "«.<■&.< + E ”’> D : 


(50) 


#=6 


#=M 


Overall energy: 

p(l7— + K— + W^) - (tf£2 + + W% 

n d( di) 8( ; 1 8£ 8^ ac' 

15, 3T. 1 d ,02;, 

“ Re d^ Kh d(^ + Re 0( ** + ** d( ^ 

1 Kfc d^gdT 1 5 (/c„ + Ke) d^/gdT v 
Re 0( d( Re dC y/g d( d£ 
p r d, TT dU, d ,. r dU. d /TT dV , 5 /Tr <9V xl 

+ #e <9(^ 0( ^ + 512 0( ^ + ^ 2 ^ + ^ 

rirr r » 

P 


, . rr r ^{7 ^5ll ^012 rr^ ^ 12 . t/^ ^022, 

+ j? e [ ac + ac <9( §c J 


1 5 , P-Le, L <96*, . 11 

D. D_ * Q/* ) 


™ pLe, dy/gdC, 
RedC ^ "" ' ' Re y/g Pr 0( d( 


(51) 
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Species continuity: 


TT dC, ,rdC. dc. 

r u -B( +pV -^ " 

1 d fiLe.dC. 1 1 iiLt'd^gdC, 

~ ~Re~dC Pr dC Re ^fg Pr d( d( 

+ w* (5 = 1,2, . . . , 11) 

(52) 

The equation of state is 


9ft Too N C$ ^ C e ^ ^ 

p - (yr> -ST Ml *’ 

U <x>M 00 *=1 e 

(53) 

where 1 is the universal gas constant, and M«, is a mixture molecular weight 
at the freestream. 

Similarly the equation of state for electrons is given by 

, Wco y G e 

p - = { vmJ p Mi T ’ 

(54) 


In equation (50), ey is the vibrational-electronic-electron energy of mixture 
per unit mass, and h v ,. is the vibrational-electronic-electron enthalpy for 
species a per unit mass. These are defined by 

ey = e„ + erf (55) 

hy,, — h v>t + hei,« 


The enthalpy of the mixture, h, is expressed by 

h^C.K 

»=i 


where h, is the enthalpy of species s including 
Supplementary relations are given by 

11 


= 1 


the heat of formation. 


(57) 


(58) 
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(59) 


11 


= 0 


«=1 


The normalized heat-transfer rate to the wall is given by 

/iLe, u dCc 


g w = U - («. + - LA -fT’ 1 dc, 


dT Vl 


(60) 


The non dimensional quantities appearing in the above equations are the 
Reynolds, Prandtl and Lewis numbers defined as: 

_ Poo^oo-^re/ 

Mr e/ 




Pr =1 Ezl 
V 


pV.Cj, 

ie - s_ 5 T 


_ pDjCpo 
Levs — _ 

A. tj 


where c p and are the frozen specific heat of the mixture and the specific 
heat of the mixture for vibrationally excited molecules, respectively. 

In the governing equations (46)-(52), the specific second-order terms are 
underlined. These terms become negligible when the radius of curvature is 
very large compared with the local boundary-layer thickness (first-order the- 
ory). Other second-order effects are present in the momentum and energy 
equations through metric tensor components which are influence yt e 


27 



* ^ it, turn affect the Christoffel symbols, and the con- 

surface curvature and which ^ related to the general 

travariant velocity components V ^ momentum equation 

base vectors a . It must be noted a surtace no^ ^ 

(49) contains a ^"“'curvature. The pre«™' » «“• Mns * ant 

gal otce ' n "' d y , In the first-order approximation, the metrics can 
across the boundary layer. significant variation across 

he estimated on the body surface smce there ., no sign. ^ ^ of 

"n’wiU he significantly reduced, 
layer equations (46)-(52) witn me . coor( Jinate system. In the 

tions [41-44] in a non-orthogon^ to J ODd order in 

viscous-shodt-layer - « retained in the 

£ ’ Ld steady Navier-Stokes equations. Also, all viscous terms in the 
normalized steady wavier h eUminate the cross derivatives 

normal momentum equation g equations are the same 

that appear. The resulting viscous- shock-layer (VbL j eq ^ ^ ^ 

to the second-order additionally applying thin layer 

face normal momentum q \ - ai^c hvDerbolic-parabolic 

approximation to the normal momentum tol ^ P „ aboKc n ature 

nature of the VSL equations can he changed to the totally par 
of the second-order boundary-layer equations 145]. 

The perfect gas boundary-, aye, equations can be «^"s du- 
ping vibrational-electronic-electron energy equation (50) and «P 

P TeSned In 

boll ^Thltegeometricqu antities are then substituted into the boundary-layer 
equations (46)- (52). 


2.4 Boundary Conditions 
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The boundary conditions at the wall axe prescribed by the usual non-slip 
conditions. The velocities at the wall are given by 

UitnC = 0) = = 0) = W((,v,C = 0) = 0 (61) 


The translational-rotational temperature at the wall may be given by ei- 
ther the prescribed wall temperature or the prescribed heat flux through the 
surface. That is: 


Ttf.n.C = O') = TJLti) 


ffi2i 


or 


0T(t,V,C = 0) _ d TJ^v) 
dC dC 

The T t = T v approximation may not be valid in the boundary layer adjacent 
to the wall, where electrons are affected by the electrical field produced by a 
charge separation (the plasma-sheath effects) [46]. Consequently, the electron 
temperature T e and the electron pressure p e at the wall must be derived using 
the plasma-sheath theory discussed in Ref. 46. However, we tentative y 
ignore the plasma-sheath effects for simplicity and assume T e = T v at the 
wall. Then the vibrational-electronic-electron temperature T v at the wall is 

™,.,C = 0) = T^) W 


or = o) _ ar,(i,i7) (6b) 

ac d( 

The body is assumed to have a fully catalytic or noncatalytic wall. Then, 
the species concentration at the wall can be expressed by 

= 0) = {for fully catalytic wall) (66) 


aC«(£,7?,C £) _ q. (far noncatalytic wall) (67) 

dC 

Usually, C tw {Z, tj) are not known a priori and therefore must be determined 
from the species mass flux balance at the wall. However, for simplicity, 
C U v) are chosen here to be equilibrium values with a specified tempera- 
ture and a pressure. For more general catalytic boundary conditions for t e 
species conservation equations, see Ref. 32. 
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The outer boundary conditions are derived from the matching between 
the viscous and the inviscid profiles at the boundary-layer outer edge. In 
first-order theory, this leads to the assumption that the inviscid flow across 
the boundary layer is constant and equal to its value at the wall. 

V{i,ihC = S) = U*{t,ih( = S) 

= M^»(=°) (68) 


V{t,r,,<; = 6) = V' d {l,ih( = 6) 

= v ed (t,i i,C = o) ( 69 ) 


In second-order theory, the matching between the boundary layer and the 
external inviscid flow must be fitted up to the second order. Monnoyer [15, 
16] deduced the following equivalent inviscid velocity distribution across the 
boundary layer assuming that the outer inviscid flow is irrotational: 


17(£,tj,< = S) = g u a u l = 0) ( 7 °) 

y(*,77,C = 6) = g 22 a 22 V td ((,v,C = 0) ( 71 ) 

The relations (70) and(71) reflect the skewing of the inviscid velocity flow field 
across the boundary layer. That is, the velocity vector 2U d is not contained 
in one plane normal to the surface as in the first-order theory Instead the 
angle between the x 1 line at the wall and y* d varies across the boun ary 
layer The rationale of these arguments can be found in Ref. 15. All other 
flow variables at the outer edge such as temperature, pressure and species 
concentrations are interpolated from the inviscid solution profiles. 


2.5 Displacement Thickness 

The displacement thickness 6 l is an important boundary-layer parameter 
in the evaluation of the extent of the viscous flow field. Assuming that the 
boundary-layer thickness <5, as well as the velocity distribution along the 
normal to the wall are known at any point on the body surface (strictly 
speaking, the boundary-layer thickness must be determined m the solution 
process by a condition of smooth merging of boundary-layer flow into the 
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outer inviscid flow), the displacement thickness is defined in such a way 
that the mass flux in the presence of the boundary layer is the same as e 
mass flux in a inviscid flow terminating at <5i instead of the wall [15, 16J. 
This definition leads to the relationship: 

[ M\p ed U ed - pU]d() = M PtJ U e<i dC) (72) 

A(v6jf* M[p ed V ed - P V}d() = M Ped V ed dC ) (73) 

where a = det(ctap) = a n a 22 ~~ ( a i 2 ) 2 - 
M is given by [19] 

M = 1 - (Ki + K 2 )C + KiK 2 C ( 74 ) 

where K\ and K 2 are the principal curvatures of the surface. Equations 
(72) and (73) are integro-differential equations for the unknown displace- 
ment thickness 6 X which must be solved iteratively after all the flow variables 
are known. The displacement effect of the boundary layer on the external 
inviscid flow can be dealt with by means of an effective thickening of the 
body, which requires the body geometry to be readjusted during the viscous- 

inviscid matching process. 

An alternative way of imposing the displacement effect is the equivalent 
source distribution proposed by Lighthill [47]. For the surface C = *1 to 
be an inviscid streamline surface, an equivalent wall outflow {p e W e )w has to 
be applied which obeys the relationship [15, 16] : 

Mp«iW' d ) w = Mp ' dUedd<:) (75) 

y/a(p ed W ed ) w = ^(t/® J Q Mp ed V ed d() (76) 

Combining equations (72), (73), (75) and (76) provides the following rela- 
tionship which no longer includes 5^: 

<w -ww. = 2^4 (v/S jf MWdV ' d ' pm) 

+ [ S M[p ed V' d - P V}dC)} (77) 

OT} JO 
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The equivalent source distribution can be calculated from equation (77) once 
the boundary layer solution is obtained. The influence of the boundary layer 
on the inviscid flow is then accounted for by a distribution of equivalent 
source distribution at the wall (77). 

A matching of the boundary-layer flow to the external inviscid flow can 
be achieved by an iterative process in which the outer flow is modified by 
the newly computed displacement effect. The resulting equivalent inviscid 
flow is then further coupled with the boundary layer by providing the new 
boundary conditions at the boundary-layer outer edge. In addition, the other 
second-order effects such as the entropy and total enthalpy gradient effects, 
which originate in the inviscid part of the flow, can be taken into account in 
the calculation by matching process. Consequently, in the future paper it is 
necessary to address the numerical methods for the boundary-layer solution, 
the computational algorithms for the inviscid Euler solution, and the viscous- 

invicid techniques. 


2.6 Thermodynamic and Transport Properties 

In the calculation of viscous flow over hypersonic vehicles, the correct eval- 
uation of the heat transfer rate to the body surface is strongly subject to the 
accuracy of the thermodynamic and transport properties which involve un- 
certainties in the high temperature regime especially in the thermal nonequi- 
librium multi-temperature environments. The nonequilibrium flow environ- 
ments surrounding these vehicles will significantly impact the aerodynamic 
and thermal loads to the vehicles. Consequently, aerothermodynamic evalu- 
ation under these circumstances requires a reasonable model for these physi- 
cal properties. Lee [29] developed an explicit formula for the nonequilibrium 
thermodynamic and transport properties based on the three-temperature 
concept (translational-rotational, vibrational, and electron-electronic excita- 
tion temperatures). To derive the nonequilibrium thermodynamic properties, 
it is assumed that 1) the rotational mode of molecules is fully equilibra e 
with the translational mode of heavy particles and that 2) the population 
densities of vibrational and electronic excitation energy levels have Boltz- 
mann distributions with a vibrational temperature T v and an electron tem- 
perature T e . In addition, an assumption is made that there is no coupling 
of energy level between these modes. The nonequilibrium transport prop- 
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-rties are evaluated by assuming that all distribution functions are close to 
Maxwellian about their inherent temperature, that is, T for eavy par i 
(molecules, atoms, and ions) and T e for electrons. 

2.6.1 Thermodynamic Properties 

By slightly modifying the formulation of Lee [29], the vibrational energy 
of the mixture e v can be expressed by 

_ _ EnLi «,« (78) 

e '- p 

where e„,, is the vibrational energy for species s given by 

rT* 

r«f 


e... = £ c^.J f 


(79) 


where c.. . is the specific heat of species s at constant volume for vibrational 
excitation, and T is the dummy variable of integration It must be no e 
that e Vt , is zero for atoms and electrons ( s - 1,2,0, /, nj. 

The electronic energy of mixture e«j is given by 


Cel = 




(80) 


where e eI> . is defined as the electronic excitation energy for species s except 
for electrons. For species s, e. 1 ,. is expressed by 

(81) 


e„, = £ = jf’ Wfl' 

JTr.f JT ”t 


where is the specific heat of species s at constant volume for electronic 
excitation. 

For electrons, 5.,,. is defined as the translational energy of electrons, i.e„ 


- _ % / $ \Tp _ ^ / & \rp 

*''•'=*'-2 { W e )7e 2 { M t )Tv 


(82) 
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The vibrational enthalpy of species s, is identical to the vibrational 

energy e v , for molecular species s. ^ _ 

The electronic enthalpy of species s, h e± ,., becomes identical to c el ,. for aU 
species except electrons. For electrons, is S lven 


he J>e = e ei , e + (^)T e -^ ( Fe )7„ 


(83) 


The enthalpy for species s, h. y for the present two-temperature model is given 

by 

T5T _'T 

(84) 


Cpej mdT + h t 


+ Z r- )iT + + Jt> 

, C c m , and c«i . axe the specific heats at constant pressure for species 
s for translation, rotational excitation, vibrational excitation and electronic 
excitation respectively, and H° t is the heat of formation. 

The frozen specific heat of the mixture c„ , the specific heat of the mixture 
for vibration ally excited molecules and the specific heat of the mixture 

for electronic excitation Z v e i axe defined by 


11 


»=i 

li 

Cpu = 53 CiCpo,* 

i=l 

11 

Cpel = 53 


(85) 

( 86 ) 
(87) 


J = 1 


Sine both the translational and rotation J modes are assumed to be fully 
excited, the specific heats for those modes therefore reduce to 129) 

, , = l JL (88) 

*** 2 M. 

for the translational mode of species s including electrons, and 


_ JR_ 
” M, 


(89) 


34 



for the rotational mode. 

The internal modes separation assumption is valid only at low tem- 
peratures. At sufficiently high temperatures, the coupling effects between 
rotation-vibration modes and rotation- vibration-electronic modes become 
important when evaluating the internal partition functions from which the 
corresponding specific heats can be derived. Balakrishnan [48] evaluated the 
partition functions for internal modes by introducing correction factors an 
generated the following curve fit formula for the vibrational and electronic 

specific beats: 


,4.186 • 10 7 , -5- r? 

Cvv,i — ( ^ H-'V* + + jfl > 


(90) 


4.186 • 10 7 ..-r -p rp . 


(91) 


^jF A'***.' 1 ^ 

where ~A, ~B, and V are constants in the above polynomial equations and are 
presented in Ref. 48. 

For tie translational mode, tie following thermodynamic relation applies 
[29]: 


For internal modes, 
where z represents r, u, or el 


S? 

Cpt, ‘ _ ^ + M. 

*pz t , ~ ^VZ % * 


(92) 

(93) 


2.6.2 Transport Properties 


In Ref. 29, tie transport properties are evaluated by extending Yos s 
formula [49], which is based on the first Chapman-Enskog approximation, 
to the three-temperature gas mixture. In Lee’s analysis [29], the controlling 
temperature in the collision integrals is assumed to be the electron tempera- 
ture T t for collisions involving electrons, and the heavy-particle translational 
temperature T otherwise. By setting T t = T v , the necessary transport prop- 
erties can be obtained for the present two-temperature approximation. 
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In the subsequent evaluations of transport properties, two modified col- 
lision integrals, Sfl’t T) md A ^(T), and the molar concentration of species 
s, 7„ are extensively used. These quantities are denned by 

(94) 

(95) 

(96) 


_ if 2M *^ r ]!/2 

A " (T)_ 3 i 7rW(10 r . + M r ) J 

A< 2) m - — [ 

A - (r)_ 5 Kwt (M. + M r y 

- = -^ = S 

7< pM, M. 


The collision integrals irlS™ and ir^ for an 11-species ^r model are 
evaluated and curve-fitted as a function of temperature in Ref. 50, and 
tabular version of the data are presented in Ref. 40. 

The viscosity of the mixture is given by 
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«Vf, 




(97) 


as 


The translational thermal conductivity x, of heavy particles is expressed 

(98) 


*1 


| ^ 10 ## 

r k S Bil M,C(m s.MT.s'j’ir.) 

where a, r is defined by 

[1 — (Tn,/m r )1 [0.45 — 2.54(7n,/Tn. r )] 

°.. = 1 + L [T+ (m. /mT iP 

The rotational thermal conductivity of mixture k t is expressed as 


7. 


«r = k S 


7. 


h Eiil T^CO + TeA^) 


(99) 


(100) 


Thus, the frozen thermal conductivity of the mixture K h for translational- 
rotational energy of heavy particles is given by 


7t h = K t + *r 


( 101 ) 
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The vibrational thermal conductivity of the mixture k v is equal to the rota- 
tional thermal conductivity k t . That is, 

k v = Hr ( 102 ) 

The V-V thermal conductivity T? v in Eq. (49), which originates exclusively 
in molecule- molecule collisions, can be written as 


K v ~~ Kv 53 


S)r=Af 


(103) 


(104) 


.iS£”,7,£?(7) + 7.£?( r .) 

The electron thermal conductivity K r is given by 

15— 7 

~ 4 k V r i 1 1-45% A^T,,) 

The thermal conductivity of electrons J? e due to collisions only between elec- 
trons is given by 

— A * ( rn \ 

(105) 


K t = K e - 




The binary diffusion coefficient of an s -r pair of heavy particles is given 


by 


kT 


” ’ fspcn 


(106) 


where p is the pressure. 

The binary diffusion coefficient between electrons and heavy particles is ex- 
pressed as 

T5„ = -J£=- < 107) 

P^(Tv) 

The effective diffusion coefficient of species s is then defined as 

jr _ ~ M «T>) ( 108 ) 

where the total molar concentration of the mixture % is defined by 

7, = ET. ( 109 > 


»=1 
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For ions, the ambipolar diffusion coefficient which has the double vdue of 
the ionic diffusion coefficient D„ is introduced due to the ambipolar diffuston 

approximation [29]; that is 

TT,=iE. (a = 6,7,8,9,10) ( U °) 

Then, the effective diffusion coefficient of electrons D, is expressed by 

7, _ glkSSt (111) 


2.7 Chemical Kinetic Model 


Under chemical and thermal nonequilibrium environments, there are cou- 
pling phenomena among the internal energy modes and between chemical 
processes and the internal modes. Therefore, the chem.cJ reaction rate co- 
efficients are no longer a function of a single temperature. Furthermore, 
several energy reactive source terms and physical properties originated in 
The enerrexchange processes appeal in the vibrational-electron, c-electron 

energy equation (34) or (50). 

2.7.1 Chemical Reaction Rate Coefficients 

Park 136 371 proposed a controlling-temperature concept for certain types 
of reactions and asLsed its validity by comparing 
with the available experimental data fo^ nitrogen andair. H gg 
using a geometrical average temperature T„ as a contro ng empera 
dissociative reactions which implicitly accounts for the vibrat.on-dissocation 
coupling. The average temperature T a is defined by 


T =T , T <1 “ ,) (q = 0.5 or 0.7) 


( 112 ) 


He also suggests that the controlling temperature for the electron impact 
ionization H T, (T. for a two-temperature mode,). All other reac *« m» 
characterized by the translational-rotational temperature T. Park recom 
mended the rate coefficients for the chemical reactions in the nonequib num 
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high temperature air be expressed as a function of the controlling tempera- 
ture, and tabulated forward reaction rate coefficients k fiP for 11 air species 
with 33 reactions [51]. The backward reaction rate coefficient k biP is related 
to the forward reaction rate coefficient through the equilibrium constant K e>p 

in the form •=- 

(113) 


k b , p 


3T e . p 


Park [52] generated the empirical formula for the equilibrium constant K e>p 
with polynomial functions of translational-rotational temperature T. 


2.7.2 Reactive Source Terms 

The effective collision frequency of electrons V* tr , which appears in the 
elastic energy exchange term between electrons and heavy particles m Eq. 
(34) or (50), is explained in detail in Ref. 29. 

The vibrational energy reactive source term w t D t represents the rate of 
change of vibrational energy of the diatomic molecules due to dissociation or 
recombination. The electronic energy reactive source term n e accounts 
for the rate of electron energy loss due to electron-impact ionization. These 
two terms are discussed in Ref. 40 and the tabulated values are presented 

there. 

2.8 Vibrational Relaxation 

The rate of change in the population of the vibrational states at low 
temperatures is described well by an equation of the Landau-Teller form [54] . 
However, at high vibrational temperatures the vibrational energy exchange 
process is primarily governed by a comparatively slow diffusion process [ J. 
A bridging function between the faster Landau- Teller relaxation rate and the 
slower diffusive rate has been proposed by Park [36, 37, 53]. The resulting 
expression for the translation-vibration (T-V) energy exchange term m q. 
(34) or (50) is 



where the bridging function, S „ is given by 

S, = 3.5 exp{—B.jT t h) 


(114) 


The quantities T. h and T Vi , h aie the translational-rotational and vibrational 
temperatures ^immediately behind the bow shock wave. The characteristic 
temperature, 6 # , for molecular species 3 axe listed in Ref. 55. 

The relaxation time of species s, T. is given as a combination of the 
Landau-Teller type relaxation time and the collision limited relaxation time 

such that 

t. = <t.>+t c . (115) 

The average relaxation time < r . > is defined as [29] 


< T, > = 


E 10 

r=l T tr 


(116) 


where r. T is the T-V relaxation time for molecular species s with collision 
partner r, and is usually determined by the semiempirical correlations of 
Millikan and White [54]. The collision limited relaxation time T ct is given by 

[36, 37] 

r = -W (H7) 

c ' c.<r v N. 

77. is the number density of species s, and c, is the average molecular speed 
of species s, expressed by 

c. = (SWT/k'M .) 112 (ii g ) 

The expression for the limiting collision cross section, <r v , is assumed to be 
^ = HT 21 (50,000/T) 2 m 2 (119) 
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3. Stagnation-Point Flow 


3.1 Basic Equations for Stagnation Flow 


If the origin of the coordinate system is at the stagnation point, then the 
boundary-layer equations can be solved by marching away from the stagna- 
tion point. One must first determine the stagnation point location and then 
needed grid. This can be done by solving the Euler equations in the Cartesian 
coordinate system. The stagnation point location is determined by interpre- 
tating the inviscid flow properties between Cartesian grid points. The origin 
of the coordinates for the Euler equations then is shifted to the obtained 
stagnation point location. At the same time, the boundary-layer coordinate 
system originates from the stagnation point is considered [56]. This technique 
has been successfully used in the axisymmetric analog method [3, 4, 57, 58] 
to compute approximately the three-dimensional boundary layers. In this 
approach, the surface streamlines which start from the stagnation point are 
calculated by using either surface velocity or pressure obtained from the 
inviscid flow solution. To ensure the efficiency and accuracy in the viscous- 
inviscid coupling, it is desirable to choose the new Euler grid after the origin 
shfting coincides with the boundary-layer grid on the body surface. Thus, 
interpolation between the Euler and the boundary-layer grid points becomes 
necessary in surface normal direction only. 

In order to start the boundary- layer solution from the three-dimensional 
stagnation point, it is necessary to obtain the stagnation-point equations 
and their solutions first. We let the origin of the surface-oriented coordinate 
system coincide with the stagnation point S on the body surface. 


That is 

s = s((= v = c = o) 


( 120 ) 


It is assumed that the stagnation line is a straight line normal to the body 
surface. In accordance with the approach of Hirschel and Kordulla [19], we 
assume further that the body surface is not too strongly curved in the vicinity 
of the stagnation point, and that the body and the flow past it, and therefore 
the flow variables, are symmetric in regard to two planes (£ = 0 and rj = 0) 
at the stagnation point. 


£= 0 , 77 = 0, C> 0: 
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dU_ 

dr] 

dT 

dt 

dp 

dt 


dV' 


= 0 


de\ 


= 0 


dt 

dT dT v _ dTv dh _ dh _ dey 

~dv = ~dt ~ ~dv ~ dt ~ dr] ~ dt dr] 
dp _ dp _ dp _ dp e _ dp e _ dC _ dC, _ ^ 
dn = dt ~ dv ~ dt ~ drj ~ dt dr] 

The definition of a three-dimensional stagnation point is given by [19] 

£ = 0, 7? = 0, <>0: 

U* = V* = 0 

Here, new dependent variables at the stagnation line (£ = 0, t? = 0, C > 0) 
are introduced 

(123) 


( 122 ) 


*«>-£ A! < c)s ^ 


Since the variables W, A', A’, p, p, ft, T, T„ and C. are functions of only 
( along the stagnation line by definition, the following set of ordinary differ- 
ential equations can be obtained by applying the conditions (121) and (122) 
to the boundary-layer equations (46)-(52): 


Continuity: 


^-momentum: 




A A 1 


(124) 


p(A 1 ) 2 + pW— = P'Mld) 2 


1 d , dA\ 11 dA 1 d^/g , V 
+ ——(p-^-) + Turf j( d( + iie 31 d< 


Re i( 


(125) 


77- momentum: 


A A 2 


p{A 2 ) 2 + pW-^r = ped{A \ d ) 2 

1 d , dA 2 ^ 1 1 dA 2 dyjg 2p 2 dA 

+ ~Redt^ dt ) + Re dt + El - 4L 


(126) 
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C* moment um: 


1 = 0 


(127) 


Vibrational-electronic-electron energy: 


+ 


J ' i?eJC <£C # 

1 1 11 


J- jL(^£l±!£ h + v V^ d^jdC. 

Redd ,=i -P»\> '* <fC -Pr p V,< d£ 




*—M 

10 


T ^-t' |5,_1 + 3 (^~)P(T ~ Tv) £ 7 ^ 


- + E 


#=6 


i=M 


(128) 


Overall energy: 


p w^ = ±±(«w ) , 1 «n‘V?<n’ 

' d( Re d( > + 

+ ££[<«. + *.)!£] + > 

726 5 c Re yjg dd dd 

_i i ^CY' —L\ i 1 1 ^,/ii/e, yd^/gdC, 

Redd^ Pr * dC + ^ 129 ) 

Species continuity: 


p\Y ^—L — 1 1 z /il/e, dC t . 

dc Redd { Pr dd } 

I 1 1 /lire, dC ' t 


Rey/g Pr dd dd + (-• 1,2,... ,11) 


(130) 


The equations (125) and (126) are obtained by differentiating the momentum 
equations (47) and (48) with respect to ( and 77, respectively. Furthermore, 
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the pressure gradient terms in the equations (125) and (126) have been re- 
placed by the convective terms at the outer edge of the boundary layer. 


Pedi-A-lj) 2 - 


922 d 2 p 9u &P 

g d(d( g d£dr) 


(131) 


, A 2 x 2 _ 9n d 2 p _ gu d 2 p 

Ped{A ed ) - g g 


(132) 


The obtained (6 + NS) equations (124)-(130) plus the equations of state 
(53) and (54) would determine the (8 + NS) unknown variables W , A , A , 
Pi Vi Pet T , T v and C,. 

Before solving the equations (124)-(130), a transformation of the normal 
coordinate may be introduced by changing the independent variable from C 
to x ' m such a way that 

X = | ( 133 ) 

where 6 is the local boundary-layer thickness. This transformation, which 
does not affect the geometrical considerations, gives the effect of fixing the 
body surface at x = 0 411(1 tlie outer houndary-layer edge at x = 1, and 
a finite number of grid points can be used over the computational domain 

0<X<1- 


3.2 Boundary Conditions for Stagnation Flow 
The usual boundary conditions at the wall are 


W(c = 0) = 0 

(134) 

A\C = 0) = >*« = 0) = 0 

(135) 

II 

o 

II 

(136) 

dT{( = 0) dT w 

(137) 

d( d( 
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(138) 


T„« = 0) = T. 


or 


dT„(C = 0) _ dT w 
d( dC 

C,(( = 0) = C tw \ (for fully catalytic wall) 


(139) 

(140) 


or 


— — = 0; (for noncatalytic wall) (141) 

dC 

The velocity W and the partial derivatives of velocities A 1 and A 2 at the 
outer boundary are 

W(C = S) = W ed (142) 


A\C = *) = A] d (C = S) = g^arrAUC = 0 ) 

A 2 (C = S) = A 2 ed (( = S) = g 22 a 22 A 2 ' d (C = 0) (143) 


where 


A\ d (0 = 


dU ed 

d( ’ 


^(0 = 


dV' d 

drj 


The outer boundary conditions for pressure, temperatures and species con- 
centration are given by the external inviscid flow solution at the stagnation 
line. 


3.3 Displacement Thickness for Stagnation Flow 

By imposing the conditions (121) and (122) into equations (72) and (73), 
the relation for the displacement thickness along the stagnation line becomes 

[* M[p,d(A\ d + A\ d ) - p(A l + A 2 )]dC 
J 0 

= [*' Mp e d(A\ d + A\ d )dd (144) 

Jo 
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Similarly, an equivalent wall outflow at the stagnation point is obtained 
by relation (77) 

= \j‘ + A h) - p( a ' + < 145 ) 

Since M is defined by equation (74), the equivalent source at the stag- 
nation point can be calculated from equation (145) once the boundary layer 
solution is obtained. 
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4. Concluding Remarks 


The second-order hypersonic three-dimensional thermo chemical nonequi- 
librium boundary-layer equations in generalized curvilinear coordinates are 
derived by means of an order-of-magnitude analysis. The two-temperature 
concept is extensively applied to describe a thermal and chemical nonequilib- 
rium flow around a general body geometry. Corresponding three-dimensional 
boundary conditions are shown along with the boundary- layer parameters 
and the invicid-viscous layer matching procedure. In addition, the expres- 
sions for the nonequilibrium thermodynamic and transport properties and 
chemical reaction rates are given. The resulting set of basic equations will be 
suitable for the development of a boundary-layer computational code which 
allows simultaneous handling of both general three-dimensional configura- 
tions and general thermochemical nonequilibrium gas models in an engi- 
neering context. By investigating a systematic analysis of three-dimensional 
flow and high- temperature nonequilibrium chemi co-physics described in the 
present paper, the future users of the three-dimensional code being developed 
will be able to have an opportunity to evaluate the accuracy of the governing 
equations that are being solved numerically (code validation). 
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Appendix A 


Metric Tensors 

By applying the definition [equation (7)], the components of the covariant 
metric tensor in the present partially orthogonal system are given by: 

< 7 n = x| + + z \ > 012 — 02i == x ( x n + ViVn + 

022 = X 7 V + yl + z\ , 013 = 031 = 023 = 032 = 0 
033 = x ( + + Z \ ~ 1 — 1) 


The components of the contravariant metric tensor can be expressed sim- 
ilarly by the relation (9): 


0 


33 


a” = 221 

s 

22 _ 

0 

011022 ~ 012 _ j 

0 


0 12 = ^ 

0 

0 13 = 0 3i = 0** = 0 32 = o 

(A -2) 
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Appendix B 


Christoffel Symbols 

The Christoffel symbols of the second kind axe given here for the present 
semiorthogonal coordinate system: 

T.1 1 r d 9u _ dguv (T3 l\ 

r " = 2 J 9 "~dT - 9n(2 ir - w }1 


„ n T .1 1 r 5^22! 

r 21 “ r i2 - 2^(522 912 ] 

(B-2) 

Til 1 r /o^5l2 dg 2 2^ 8g 2 2, 

r 22 2g [522( 2 ^ ^ ) 012 ^ 1 

(B-3) 

_2 1 r in dgv & 9 ii\ 8gu 1 

Tn ~ 2g^ 2 d( dr] ^ 913 d( ^ 

(B — 4 ) 

_2 r >2 1 r ® 9 22 ^011 1 

r 2i - r i2 - 2 g [ 0 11 di 912 drj^ 

(B — 5) 

t,2 1 r dg 22 ,089x2 8922 w 

r “-2 g l9u ar, - 911(2 ft, - * 11 

(B — 6) 


(B-7) 

T .1 r>l 1 r 0011 . 00 12 1 

1^13 “ r 31 - 2g [922 ^ 01*^1 

(B — 8) 

nl Til If 0012 - 00 22 1 

r 23 - r 3, 2 - 2^1022 ^ 9X2 Ql > J 

(B - 9) 

n 2 ti2 1 r 0012 _ 00Hl 

r !3 - r 31 “ 2^(011 Q(- 9 X 2 ^ ] 

(B - 10) 

_2 n2 0022 0012 1 

r 2 3-r 3 2 - 2y [0ii d( - 512 ^ 

(B — 11) 

= r|,i = o; r‘ 3 = o ; (* = 1 , 2 , 3 ) 

(B - 12) 


where g = 011^22 - (^n) 2 - 


49 



The co variant curvature tensor b a p of the surface is related to the Christof- 
fel symbols of the second kind by 

tw = i%>; K/S = i,2) (B-is) 

where the dot superscript *•’ denotes a quantity on the body surface. From 
equation (B-7), the following expressions for b a p result: 

, _ 1 dan 

bu ~ ~ 2 d( 

, 1 $<*22 

“ '2 d( 

- -rw (B - 14) 
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Appendix C 


Principal Curvatures of the Surface 

The principal curvatures K x and K 2 [equation (45)] play an important role 
in second-order boundary-layer theory since they indicate if the local radius 
of curvature of the surface is not negligible compared to the boundary-layer 
thickness. K x and K 2 are determined from the following relations [19]: 


K x + K 2 = 


On&22 ~t" 022^11 ~ 2a i2 bi 2 


Oll®22 ~ ( a 12) 2 

(C-2) 

The Principal curvature radii R x and R 2 axe given by 


(C- 1) 


„ „ b u b 22 - (bi 2 ) 2 

Jii • ii 2 — rz 

O11O22 — v° 12 / 


Rl K x 


R 2 = 


~K, 


(C~ 3) 
(C~ 4) 


The directions of the principal curvatures of the surface Ai and X 2 can be 
obtained as the solution of a quadratic equation [19]: 


A 


1,2 - 


± 


&22 a n ~ bn°22 

2(6nai2 — &i2fln) 

y/(biia22 — 522 Qh) 2 — 4 (bnai 2 — bi2Qn)(bi2 ffl 22 ~ 622^12) 
2(bnai2 — i>i2 a n) 


(c- 5 ) 
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Fig. 1 Orientation of coordinate systems. 

(y i '(x,y,z): Cartesian coordinate system; 

x ! (!,r)i£): Curvilinear surface-oriented coordinate system) 
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